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We present a model for biological aging that considers the number of individuals whose (inherited) 
genetic charge determines the maximum age for death: each individual may die before that age due 
to some external factor, but never after that limit. The genetic charge of the offspring is inherited 
from the parent with some mutations, described by a transition matrix. The model can describe 
different strategies of reproduction and it is exactly soluble. We applied our method to the bit-string 
model for aging and the results are in perfect agreement with numerical simulations. 
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C/5 ■ I. INTRODUCTION 

m ! 

Aging is an extremely complex biological phenomenon of immense importance and interest. Recent progress in 
studying aging (or senescence) has pointed out the importance of both genetic and environmental components ■ 
Although no single theory fully explains all aspects of the aging phenomena, the evolutionary and free radical theories, 
in particular, are supported by significant observational and experimental evidence. Evolutionary explanations of aging 
fall into two classes. First, according to the optimality theory, organisms might have evolved the optimal life history, 
' in which survival and fertility late in life are sacrificed for the sake of early reproduction. Second, the life history 
might be depressed below this optimal compromise by the influence of deleterious mutations; since selection against 
late-acting mutations is weaker, deleterious mutations will accumulate, i.e., impose a greater load late in life. 

Besides experiments with flies and data from human populations, computer simulations are a widely used tool to 
(3JQ 1 study aging . The recently introduced bit-string model , based on the mutation accumulation hypothesis jl]-[| , 
is able to reproduce the exponential increase of the mortality with age || - this behavior is known as the Gompertz's 
law. Some applications of this model have been reviewed for both asexual and sexual versions, with emphasis on 
i' computer simulations. However, only a few analytical results on this model are available. Namely the exact results 
for the survival rates Jl0|,[llJ and dynamical aspects (lj] of the catastrophic senescence of the Pacific salmon and the 
description of stationary states through Leslie matrices jl3| . In this work, we present analytical results for a general 
asexual model for different strategies of reproduction. Even for the special cases where analytical results for the bit 
string model are available, our techniques are more easily and efficiently implemented. This paper is organized as 
follows: in the next section, we present a general formalism to biological aging. In section 3, we review the bit-string 
model and build the mutation matrix for it. In section 4 we present comparisons with computer simulations. Finally, 
d ' we present our conclusions. 
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II. THE MODEL 



Consider a population with N(t) individuals living in an environment with finite resources at time t, and with 
diverse genetic charges which determine different limit ages m, beyond which the individuals cannot survive. Suppose 
the natural resources allow a maximum population N max , to be considered in a logistic (Verhulst) factor. We note 
x(a, m, t) as the relative number of individuals with age a and programmed death age m at time t, that is, 

*(a,m,t) = %M. (1 ) 

Also, consider the initial and final reproduction ages as R and Rf. The evolution of the population in a discrete time 
is described by 

x(a + 1, m,t + 1) = [1 — x(t)] x(a, m, t) for 1 < a < m — 1 

x(a+ l,m,t+ 1) = fora>(m-l) (2) 
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x(l, m,t + 1) = b[l - x(t)] ^2 Amm' X! X ( a ' TO '' *) ' 

m' a—R 

where 

m— 1 

x(t) = ^2 x{a,m,t) , (3) 

m a— 1 

6 is the number of offspring in each reproduction and A mm i is the birth matrix, that gives the probability of an 
m'-parent having an m-offspring. When no mutation is allowed, A mm i = S mm '. We also considered a — 1 to be the 
first year of life, so that a > 1. Note that the maximum age a reached by an individual is m — 1. 

Eqs.(||) describe the evolution of populations whose transmitted genetic charge contains the information about the 
maximum age of death; the way this genetic charge may change from parent to offspring through mutations as well 
as how the maximum age of death is distributed among the population depends on further details, typical of each 
species or theoretical model. The solution to equations (Q) strongly depends on the form of the birth matrix A. 

Hence, to completely investigate the evolution problem given by equations (||), we must first choose a transition 
matrix A. For that we consider in the next sections the bit-string model for biological aging, that seems to have 
grasped some of the interesting features of the age structure of populations. 

On the other hand, when this transition matrix satisfies some general features, the present model is exactly soluble, 
and before applying it to the bit-string model, we present this exact stationary solution. 

From eqs.(2j) one obtains 

x(a, to) = (1 — x) a 1 x(l, to) for 1 < a < m , (4) 

where we explicitly assume the stationary solution and from now on will not write the time dependence. Our methods 
thus do not deal with mutational meltdown where the population is doomed to extinction. 

We consider the stationary solution and the mutation matrix A mm i to fulfill the following assumptions: 

1. x(l, to — ► oo) = 0, what is observed in living populations since offspring with unlimited expected life length are 
not possible (in the bit-string model this condition is equivalent to x(l,B + T) = 0, where B is the bit-string 
length and T the lethal number of accumulated diseases); 

2. A mm i is a triangular matrix such that A mm i = for to > to', that is, parents cannot give birth to offspring with 
larger life expectancy (which corresponds to only bad mutation in the bit-string model) . This is not biologically 
unrealistic for well adapted populations - advantageous mutations are expected to be extremely rare due to the 
large times required for noticeable species evolution; 

3. A mm 7^ 0, that is the probability that the parent gives birth to offspring with the same expected life length is 
different from zero. This condition is also expected in biological populations and 

4. A mm < A m i m i if m > to', the probability that a parent gives birth to offspring with the same expected life 
length to decreases with to. In other words, the larger the parent expected life length is, the larger the probability 
that a difference in the genetic charge of the offspring effectively reduces their expected life length ( the better 
the genetic code, the larger the number of events that can spoil it). 

The solutions to eqs.(||) depend on the reproduction features of the populations. We consider two cases: the first, 
when the individuals reproduce only once, that is Rf = R, and the second case, when Rf > R and the individuals 
reproduce every year after the reproduction age, until they die. 



A. Case R f —R 

For semelparous populations (that reproduce only at age R), the relative number x(l, to) of m-babies can be written 

as 

x(l, m) = 6(1 - x) A mm >x(R, to') . (5) 

m' 

We first observe that as x(a, to) = for a > to the individuals with programmed death age m below R + 1 do not 
reproduce. Nevertheless these individuals keep being born from parents with to > R + 1, due to mutations. 
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In a population we assume a maximum life expectancy m = v (in the bit-string model discussed in the next section, 
v can be taken as B + T) and we assume x(l, m) — for m > v. Now we can write eq.(^|) for m = v: 

x{l,v)=b{l-x) R A vv x{l,v) , (6) 

since we are assuming a triangular matrix and x(l, m) = for m > v. Hence, either x(l, v) = or 6 (1 — x) A vlJ = 1. 
Now, the expression for m — v — 1 reads 

x{l,v-l)=b{l-x) R [A v _ Xv _ x x{l,v-l)+A v _i v x{l,v)] . (7) 
If x(l, v) 7^ 0, then 6(1 — x) R = A~l and the above equation reduces to 

i - - i) = y) , (8) 

which is not possible because x(l,m) > for all m and the left hand side is negative as a consequence of 
A m m < A m _i,„_i (assumption 4). Hence x(l,v) = and the solutions to eq.(Q) are either x(l,i> — 1) = or 
6(1 — x) R A„^x v ^i = 1. This situation repeats on and on up to m = R + 1, such that x(l,m) = for m > R + 1. 
The equation for m = R + 1 leaves us with two possibilities: either x(l, R + 1) = or 

6(1 - x) R A R+1 R+1 = 1 . (9) 

The first solution reduces to no population at all, since for m < R + 1 individuals cannot reproduce. On the other 
hand, the expression for x(l, m) when m < R+ 1 no longer presents the diagonal term and hence does not give place 
to non-positive solutions, that is, 

x(l,m)= f mR+1 X (1,R+1) , for to < -R+1. (10) 
Ar+i r+1 

Using eqs.(^) and (||), we can obtain x(l, R + 1): 

x(l,fl + l)= B+1 x ^13±l . (11) 

From cq.(|9|) the total population may be easily obtained: 

1 x l/R 



^ = 1- 71 . ( 12 ) 

and a critical number of offspring per reproduction 6 C may be defined from eq. ( |l2| ) : 

bc = ^ • (13) 

Ar+i r+1 

If 6 < 6 C the only solution is x = 0. The population age structure is given by the relative number x(a) of individuals 
at age a: 

x(a) = (l-x) a - 1 ^x(l,m) . (14) 



The catastrophic senescence effect [ |lO|JTT| is clearly present, since eq.([14|) implies that x(a) — for a > R since 
x(l,m) = for m > R + 1. Hence the observable quantities x and x(a) of the stationary solution are explicitly 
calculated by the transition matrix A and the parameters 6 and R. 

The birth matrix A mm / is a central point of the model. As stated before, it describes the probability of an m'- 
parent to give birth to an m-offspring, and hence it is closely related to the probability of mutations during the birth 
process. However, assuming that the birth matrix is triangular and fulfills the basic assumptions 2, 3 and 4 for a 
single mutation, then it is triangular and fulfills the same assumptions also for M multiple mutations, since in this 
later case the resulting birth matrix is the product of M single mutation birth matrices. In other words, if the above 
results apply to a model with a single mutation per birth, they also apply to M multiple mutations, provided that 
the birth matrix is taken as the product of M single mutation matrices. 
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A special case happens when no mutation at all is considered. In this case A mm > = S mm / does not fulfill the basic 
assumptions but the problem is still exactly soluble. In this case, 



XM=0 = 1 -[b) (15) 

and the age structure of the population is determined by the values assumed for x(l,m), that must only satisfy 

R+1 1 — (1 — x)" 1 ^ 1 

x — x(l,m) . (16) 



B. Case R f > R 

It is interesting to investigate what happens when the individuals may reproduce more than once, that is, for a > R. 
We will first consider that an individual reproduces from age R until its death (iteroparous population). In this case, 
birth is described by the equation 

ro'-l 

x{l,m)=b(l-x)^^mm' ^2(l-x) a - 1 x(l,m') . (17) 

rn' a—R 

Summing over a we have 

x(l,m) = b(l-x) R J2 A *m m 'X(hm') (18) 

m'=R+l 

where 

1 - (1 - x) m '~ R 

A-mm' = A mm ' • (19) 

It is straightforward to verify that since x < 1 the renormalized birth matrix A* does not necessarily conserve the 
required properties for a catastrophic senescence solution, that is, A* nm could increase with m. However in order to 
have x ^ there must exist some me such that x(l 7 me) ^ and x(a, m) = for m > me. If me > R + 1, that is, the 
senescence is not catastrophic, then the following must hold: 

A mt -i me -i [l-{l-xT^- R ] < i ^ 
Ami mt [1 ~ (1 ~ x) mi! 

so that the equation for x(1,to^ — 1) has a non null solution, in analogy to eq.(|8[). The ratio between the matrix 
terms is larger than one, and so the second ratio in the left hand side of the above equation must be small enough to 
compensate. Given an me > R + 1 that satisfies eq. (|20|) , the total population is obtained from the solution of 

(1 - x) R - (1 - x) me 1 , x 

* = Va ' (21) 

X UJ~L m ^ frig 

Observe that when me = the above equation reduces to the result obtained in the previous section with Rf = R, 

as it should. 

It is possible to obtain the maximum me for a given birth matrix A, that is, the maximum death age of an iteroparous 
population. It is useful then to define a maximum death age limit /i: 

/i = maxjmj} . (22) 

where all me satisfy eqs.(^o|) and (pl|). When starting from an initial population containing death ages larger than 
H, the solution converges to me = /x. However, depending on the structure of the birth matrix, it can happen that 
eq.(|2(]) never holds for me > R + 1. In this case, the system organizes itself as a semelparous population and all 
results of the previous section apply. Also, when the population initially contains only individuals with death age 
m < fi, then me < /i, because the birth matrix A is triangular. 
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An intermediary case, when R + 1 < Rf < /i, may also be considered. In this case the solution is rri£ = Rf + 1 for 
initial conditions where to > Rf + 1 are present. 

The solution to the problem is complete with the values of x(l,m) for R + 1 < m < mi, that are obtained by 
finding the solution to the eigenvalue equation 

1 m ' 

X(1 ' m) = l \l-(l-x)^-R] £ A mm > [1 — (1 — x) m '~ R ]x(l, m') , (23) 

for R + 1 < to < mi, that satisfies the normalization condition given by 

m — 1 

•^EE^™) ( 24 ) 

m a—1 

where xla, m) = (1 — x) a ~ 1 x(l, m). 

Here we can also discuss the case where mutations are absent. As before, A mm i = S mm >. The eigenvalue equation 
reduces to 

x(l, m) = b± - K - —x(l, to) for to > R. (25) 

x 

That means that for every to* > R there is a different solution, where the population consists only of m*-individuals 
and the total population x satisfies the above equation for to — to* . The final state depends on the initial conditions 
and the system always converges to a solution with mi given by the largest to present in the initial population. 
In the next sections we apply our approach to the bit-string model. 



III. THE BIT-STRING MODEL 

The bit-string model for biological aging consists in a population of N(t) individuals at time t, each one 
represented by a bit-string of B bits and subject to aging, reproduction and death. The bit Si (i = 1, . . . , B) of the 
string contains the genetic information of the programmed health status of the associated individuals at age i: when 
Si = 1 a genetically programmed disease is acting for all a > i. Each individual in the population may survive up 
to T of these diseases, that is, it cannot survive longer than age to, at which the sum of bits from zero to to is T. 
Individuals cannot reproduce before age R, when it gives birth to b offspring either once (semelparous) or each year 
(iteroparous) until death. An offspring has the same string as the parent except for M bits that are randomly changed, 
to simulate genetic inheritance and mutation. The environment limitations are taken into account through a Verhulst 
factor, that reduces the number of individuals by (1 — N/N max ), where N max is the maximum allowed population 
size. From now on we shall discuss in terms of the relative population quantities x — N/N max , x(a) = N(a)/N max 
and x(a,m) = N(a, m)/N max . 

Exact results have been available only for semelparous populations with T — 1 and deleterious mutations [flo|| . 
The analytical approach to aging that we present in eqs.(Q) may be adapted to reproduce the simulation model. To 
account for the finite length of computer bit-strings, we also take x(a, to, t) = for a > B whereas to > T, since 
in the bit-string model no individual may have its programmed death age below T . Moreover, to < (B + T) where 
to = B + k is associated to individuals with only k < T bits equal to 1. 

To describe the bit-string model for biological aging, we note that the probability of finding an individual with 
programmed death age to is equivalent to the probability of finding the T th genetically programmed disease acting 
on age to, or to find the T th bit equal to unity in the m th position of the bit-string associated to the individual. Also, 
the birth matrix A should reflect the change in the death age implied by the flipping of randomly chosen bits when 
building an offspring bit-string. In a recent work, Ito |T§| ] obtained an analytical approach to the bit-string model 
by considering the relative number n(S,a,t) of different bit-strings S = (Si,^, . . . ,Sb) present in the population. 
The present model reproduces the evolution equations of Ito [jl8f by summing over all bit-strings with the same death 
age to, that is x(a,m,t) — J2{§} n ( a ' S,t)5(£ — to) where £ is the position (locus) of the T th inherited disease in the 

bit-string S. Due to the form of the birth matrices of the bit-string model, the general calculations performed in the 
previous sections still apply and exact results may be produced. 
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A. The mutation matrix 



To write the birth matrix A m m < for the bit-string model, we first consider the one-mutation matrix F, when at 
most one bad mutation happens at birth in the bit-string model: a random site is chosen and set to one, regardless 
its previous state. The mutation matrix is triangular and reads 



for m < B 

TP, 

± ran 
TP 

± ran 
TP 

± ran 

for m! = B + 1 

Fran 

^ ran 
TP, 

± ran. 

for B + 2 < m! < B + T 

TP 

± ra 
TP 

1 m 

F m 
Frn 



T {ml -T)\ (to- 1)! 
B{m-T)\ (m'-l)! 
B — m' T 
B + B 



for m < m' — 1 



for to = m 



for ml + 1 < to < B + T 

T (to - 1)! (B -T+l)! 



B (m-T)\B\ 
T - 1 



for T < to < B 



B 



= 

= 
m! 



T 



B 

B + T - to' 



for to = B + 1 
for m > B + 1 

for m < to' — 2 

for to = m' — 1 



B 



= 



for 



for to = to' 
1 < to < B + T 



The matrix elements are obtained considering that the probabilities of finding a bit one in different ages before age 
to are not correlated. Now, the string length B and the maximum number of diseases T of the bit-string model have 
been already taken into account (actually B is necessary only due to the finite limits of computers. 

In the case that more than one mutation may happen at birth, say M mutations, the mutation matrix is taken to 
the power M, that is A = F . The birth matrix is triangular and considering Rf — R populations, it fulfills the 
required conditions for a catastrophic senescence for any M > 0. 



IV. RESULTS 



Considering the bit-string model, there are three different ways of obtaining results: analytically as presented in 
the previous sections, numerically iterated solutions to eqs.(^) and numerical simulations. In what follows we will 
discuss these three forms, for each case we have considered. 

We first present the results for Rf = R = 11 (age at reproduction of the Pacific salmon, in years). The birth matrix 
A is obtained from the one mutation matrix F, that is 

A = F M for M > , (26) 

and A mm i = 8 mm i for M = 0. In fig. 1 we present a snapshot of the time evolution of x for M = 1 and T = 1, 
considering different initial conditions: a) fixed point solution as given by eq.(^2|) and b) all individuals with age 1 
and to = 32, i.e, free of mutations. The iteration of eq.(||), starting from the solution of eq.([l^) confirms the existence 
of a fixed point of the dynamics where the system remains for all times. However, for any values of M and T, if b > b c 
and Rf — R we also found oscillatory states with period 11 (= R). These oscillations may vary in amplitude and 
phase, depending on the initial condition. They have been previously observed fl2] , |l4| and are due to the non-overlap 
of generations in semelparous populations |l2] ]. The effect of increasing T, the limit number of diseases that kill the 
individual, is to increase the population. On the other hand, increasing M decreases the total population, because 
offspring with low to are more frequently generated in each birth process. Simulations show that the fixed point 
solution is very unstable: any perturbation drives the system slightly away from the fixed point, after a few steps the 
population is driven to an oscillatory regime. Another point worth noting is the extremely slow convergence of the 
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age distribution. When the population size reaches the equilibrium (after less than 2000 steps) the population with 
age one year after the age at reproduction has not vanished. An exponential fitting (x(a) = Aexp(crf)) to this age 
gives a small valued exponent a = —0.004 which guarantees that the population older than age at reproduction will 
eventually vanish. 

In fig. 2 we present the survival rates, defined as 

- (^) • • * i ■ 

The catastrophic senescence effect is clearly present in this figure. The excellent agreement between computer simu- 
lations and the analytical results can be also seen. 

We consider now R < Rf, for the bit-string model. Considering B — 64 and reproduction every year until death, 
our results do not show any oscillatory behavior: regardless the initial conditions the system converges to a fixed 
point solution. In fig. 3 we present a semi-log plot of the mortality rate at age a defined as [fl5| 

for M = 1,T = 1,6 = 0.1. This normalization has been proved to eliminate the Verhulst factor influence jyj. Using 
this normalized mortality curve we can see the Gompertz's region (from ages 10 to 25). A clear change in the behavior 
at the minimum age of reproduction R - where we expect that mutation accumulation effects are not relevant - can 
be seen in this figure. Another deviation occurs also at older ages. The existence of several Gompertz's regions was 
proposed by Gompertz |i5[ and it was already studied in the bit-string framework jl6| . To obtain the analytical results 
for the asymptotic solution we must first obtain the maximum death age of the population. For these parameters, 
the maximum death age is fi — 33. Starting from an initial condition where m > 33 are present, the population 
stabilizes with m,( = fi = 33 (maximum age=m^ — 1). Solving eq.(|2l]) for the adequate birth matrices A, we can find 
the asymptotic total population. Again we present computer simulations for comparison. After an enough number of 
time steps, individuals with advanced ages (a > 33) tend to disappear from the population. Our results were obtained 
from longer series (3 x 10 e 'time steps) than than the ones from intensively parallel simulations by Meisgen |l7j (800,000 
time steps). 

It can be argued about the sizes of the simulated systems. This point is one of the important advantages of the 
present treatment. We adopted, in computer simulations, population sizes around 300000 individuals (considering 
B = 32 for comparisons, which is the most used value in computer simulations). In a Pentium 150MHz (32Mb RAM) 
running Linux and using a very optimized code, 1025 sec are needed to simulate 20000 time steps. However, solving 
the equations for this same time interval, using the present approach only just 20 sec are needed, that is, in our 
strategy is three orders of magnitude faster than a simulation (for 300000 individuals)! Also, it is worth remarking 
that no population finite size effects are present in this treatment. (The finite size effects related to the bit string 
length are exactly taken into account). Moreover the storage needed to simulate it is basically 2 * B x B x 8 bytes 
(considering double precision float point numbers) whereas a N = 300000/8 bytes were required to the computer 
simulations (using bitwise operations). However, we have to emphasize that in the current approach the computer 
time increases as B 2 whereas in computer simulations the CPU time increases linearly with B. 

In summary, the analytical calculations agree perfectly well with both the stationary solutions found by numerically 
iterating eqs.(g) and with simulations. For R — Rf, oscillatory stable solutions may also be found when solving 
iteratively the evolution equations. Finally, we observe that the computer time required for numerically iterated 
solutions for the analytical evolution equations is very short due to the size of the transition matrices (B x B), in 
comparison to previous models where the birth matrices are 2 B x 2 B p8|| ■ This CPU time required in the present 
model is still orders of magnitude smaller than that required by computer simulations of systems comparable to real 
sizes, with the traditional B = 32 age intervals. 

V. CONCLUSIONS 

We have presented a theoretical model to describe how inherited genotypes may determine the age structure of a 
population. This model takes into account the inherited maximum lifespan, which is passed from parent to offspring 
according to a birth matrix that may contemplate the possibility of mutations. When this birth matrix is triangular, 
the stationary solution is obtained analytically for semel- and iteroparous population. Considering that reproduction 
occurs for R < a < Rf, the expected catastrophic senescence is obtained for Rf — R, but it can also occur for Rf — > oo, 
depending on the structure of the birth matrix. In general cases, successive iterations (not simulations) always reach 
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asymptotic solutions to the evolution equations that does not require large computational resources. Applying this 
technique to the bit-string model, we found that iteroparous population always stabilize at the theoretically predicted 
stationary solution, but semelparous ones may also present stable oscillatory behavior with period R which amplitudes 
depend on the initial conditions. The agreement with simulation results are noticeable. This approach is at least three 
orders of magnitude faster than computer simulations and still more memory saving. It is worth to remember that 
computer simulations have been used up to now, as the most important tool to aging studies. Therefore, we believe 
that our results can be useful for further studies on biological aging. A natural extension is to consider a continuous 
time limit, where the time scale is given by R, the biological relevant quantity. This extension is now in progress and 
will be presented in due time. 
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VII. FIGURE CAPTIONS 



Figure 1. A snapshot of the time evolution of the population for two different initial conditions: dotted line refer 
to the fixed point solution of eq.([l^) and solid lines to all individuals initially free of mutations. The parameters used 
in both conditions were: Rf = R = 11, b = 10, M = T = 1. From the top to bottom the solid lines mean the total 
population, population at ages 1,5 and 12 years. We can note a displacement of the peaks at each age curve. Since 
the total population is the sum of all ages, these peaks present a smoother aspect. Population vanishes for ages above 
R = 11. However, the decay of these curves for the oscillatory solutions is extremely slow. 

Figure 2. Survival rates S(a), in units of S(l), for Rf = R = 11, b = 10, M = T = 1. Catastrophic senescence 
is clearly seen in this plot. Solid lines are the analytical results and the full circles are the results from computer 
simulations. 



Figure 3. Normalized mortality rate q(a) for R = 8, Rf = 64, M = T = 1, and b = 0.1. The full line refers to the 
theoretical predictions after 3 million time steps (this is the longest time series obtained for the model). Computer 
simulation results after 10 4 steps are represented by circles. The population beyond age 33 will eventually disappear 
(after one million years, according to the theoretical predictions). 
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